function C = C_ML_Lasso( y, S_init,lambda)

y = double(y);
S = S_init;
[height, width, wavelength] = size(y);
[~,k] = size (S_init); % k is number of pure components K是纯组分的个数
C = zeros(height, width,k);
%e = nan(height, width, wavelength);
mask = nan(height, width, wavelength);
mask(~isnan(y)) = 1;
% Initialize e as y-A(S\circ I)C  初始化
%e = zeros(height,width,wavelength);

% Initialize C as the ML estimate, with lasso as regularization
% 初始化C作为ML估计，套索作为正则化
C_prev = zeros(k,1);
sampled_sum = 0;%sample为样本
%squared_error_sum = 0;
for i = 1: height
    for j = 1:width
        non_zero_locs = find(mask(i,j,:)==1); %find函数是返回向量或者矩阵中不为0的元素的位置的索引
        non_zero_count = length(non_zero_locs);
        sampled_sum = sampled_sum + non_zero_count;
        S_temp = zeros(non_zero_count,k);
        y_temp = zeros(non_zero_count,1);
        for l = 1:non_zero_count
            S_temp(l,:) = S(non_zero_locs(l),:);
            y_temp(l) = y(i,j,non_zero_locs(l));
        end
        % Note that Lasso requires at least two rows for y_temp 注意，Lasso要求y_temp至少有两行
        if non_zero_count >= 2
            %C_temp = lasso(S_temp,y_temp);
%             C_temp = lasso(S_temp,y_temp,'lambda',lambda);
            C_temp =lasso(S_temp,y_temp,'lambda',lambda);
            %lambda为正则化系数，指定为逗号分隔的一对，由'Lambda'和一个非负值的向量组成
            %S_temp为预测数据，指定为数值矩阵。每一行代表一个观察结果，每列代表一个预测变量
            %y_temp为响应数据，指定为数值向量。其长度为n，其中n为S_temp的行数。响应y_temp(i)对应于S_temp的第i行
            %C_temp为拟合系数，作为数值矩阵返回。是一个p×L矩阵，其中p是S_temp中预测器(列)的数量，L是Lambda值的数量
            %lasso是用于正规化最小二乘回归分析的算法.在matlab中是一个可直接调用的函数.
            %lasso函数返回预测器数据X和响应y的线性模型的拟合最小二乘回归系数
            %C_temp = S_temp\y_temp;
        else
            % Take the C of previous pixel as a rough initialization 取前一个像素的C粗略初始化
            C_temp = C_prev;
        end
        C_prev = C_temp;
        %e_temp = y_temp - S_temp*C_temp;
        for l = 1:k
            C(i,j,l) = C_temp(l);
        end
        %for l = 1:non_zero_count
        %    e(i,j,non_zero_locs(l)) = e_temp(l);
        %end
        %squared_error_sum = squared_error_sum + sum(e_temp.^2);
    end
end


%Display the initial estimation of C 显示C的初始估计
%display_results_fungi;

end